Frontiers in Systems Biology

Mechanistic ordinary differential equation models are widely used in systems biology to represent biochemical networks, population dynamics, cell-state transitions, and other biological processes; however, their predictive value depends critically on accurate parameter estimation from noisy and often sparse experimental data. In this tutorial, we present the Weak-form Estimation of Nonlinear Dynamics (WENDy) method as a forward-solver-free approach that reformulates parameter estimation as a covariance-corrected weak-form regression problem by integrating the model equations against compactly supported test functions. We present the background on the methodology through the lens of the familiar logistic equation, and we demonstrate applications of the method on real experimental data through two systems biology examples: a glycolytic oscillator with relatively dense time-course data and a sparse epithelial-mesenchymal cellstate transition model with multiple experimental replicates. Ultimately, using WENDy, we estimate interpretable biological parameters with uncertainty for systems with noisy and sometimes sparse available experimental data.

Heterogeneous biomedical knowledge networks (hetnets) integrate disparate data types, drugs, genes, diseases, and pathways, across independent sources; Hetionet (https://het.io) is a widely used example. A standard approach for assessing connectivity significance is XSwap, which permutes the hetnet P times and fits a gamma-hurdle null model to the degree-weighted path count (DWPC), pooling permuted values across pairs with matching source and target degrees to increase the effective sample size. This permutation approach has been highly successful in practice, but it faces four practical constraints in large graphs: (1) a finite resolution for the smallest reportable p-values, (2) computational cost that grows prohibitive at path lengths L [≥] 4 or 5, (3) a variance model (Var {propto} {micro}2) that departs from the configuration-model form (1 +{kappa} ){micro}, and (4) O(P 10m L) runtime. To complement this approach, we present HetNetEX (Heterogeneous Network EXact inference), which computes the null DWPC distribution analytically from degree sequences using the configuration model in O(Ln) time. In simulations at P = 200 across L = 1-4, HetNetEX achieves Spearman{rho} > 0.96 concordance with XSwap rankings while being >10,000x faster and providing analytical p-values without a resolution ceiling. High-degree pairs show larger XSwap sampling error than low-degree pairs, reflecting the finite-sample nature of permutation that analytical computation avoids.

Approximate Bayesian computation sequential Monte Carlo (ABC-SMC) propagates its particles with a perturbation kernel, and with the standard Normal kernel it degrades sharply as the parameter dimension grows, a failure usually attributed to dimension itself. We show instead that it is governed by the quality of the summary statistics, with dimension entering only through a separate and milder mechanism, and that the two must act together for the Normal kernel to break. The first ingredient is covariance overinflation: the kernel covariance, estimated from the particle cloud, overshoots the true posterior covariance by a factor set by information loss in the summary statistics. We derive this overscaling factor in closed form for a Gaussian model with sufficient statistics and show that it stays modest at any dimension, shrinking toward its baseline value as the tolerance tightens; the extreme values seen in practice (of order 103) are a signature of insufficient summaries, not of dimension. The second ingredient is perturbation overconcentration: the normalised Normal step size concentrates around one as the dimension grows, so every proposal overshoots by the same factor. Either ingredient alone is harmless; only their combination breaks the Normal kernel. A Cauchy kernel (multivariate t with one degree of freedom) removes the concentration, keeping a positive acceptance rate under arbitrary overscaling at a bounded worst-case cost of 1.87x in expected squared jump distance. In a Metropolis-Hastings framework we derive closed-form acceptance rates for both kernels that illustrate the advantage of the Cauchy kernel in this limit. A series of full ABC-SMC computational experiments on five problems at d = 12, including a hierarchical gene-expression model, show the Cauchy reducing the sliced Wasserstein distance to the reference posterior by factors of up to 50 with the same simulation budget. Since the summary statistics are commonly insufficient for the models that require ABC, overinflation is structural and the Cauchy perturbation kernel is the right default for problems in higher dimensions.

The osteocyte network in bone is believed to play an important role for how bone tissues sense and respond to mechanical stimulation. Yet, bone adaptation to mechanical loads is often conceptualised as a simple response to mechanical stimuli, such as Wolffs law, which is based on mechanical variables only and takes no account of the cellular basis of mechanosensation. Wolffs law presumes the existence of a reference mechanical stimulus, the mechanical setpoint, above which bone is consolidated, and under which bone is removed. In this paper, we develop a theory of bone tissue sensing and adaptation based on osteocytes to provide new understanding of the role played by osteocyte signals in mechanical adaptation. In this theory, the mechanical setpoint of Frosts mechanostat is explicitly embodied as osteocyte properties involved in mechanotransduction. The mechanical setpoint is allowed to adapt due to the replacement of osteocytes during remodelling, making the setpoint space and time dependent. We propose a mathematical model to implement this new theory of bone adapation and present numerical simulations of this model to explore how mechanobiological response curves (effective Wolffs laws) are modulated by setpoint adaptation during remodelling. By accounting for varying osteocyte populations within bone tissue, we explore bone adaptation under osteocyte disruptions, which is particularly relevant to age-related bone loss. Our model suggests that biological disruptions of remodelling balance cannot always be compensated by mechanical feedback, and that setpoint adaptation during remodelling may have significant observable consequences, such as hysteresis in bone response signatures that resemble lazy zones.

We present the BCCWJ-Brain dataset, a multi-modal neuroimaging resource comprising functional magnetic resonance imaging (fMRI), magnetoencephalography (MEG), and electroencephalography (EEG) data recorded from native Japanese speakers reading newspaper articles from the Balanced Corpus of Contemporary Written Japanese (BCCWJ). Neural data were collected from 112 participants (36 fMRI, 35 MEG, and 41 EEG) as they read twenty newspaper articles presented in a Rapid Serial Visual Presentation (RSVP) paradigm. By providing three complementary neuroimaging modalities collected under identical naturalistic reading stimuli, this dataset provides a cognitive benchmark for computational models such as large language models. The dataset is publicly available on the OpenNeuro platform, offering a valuable resource for neuroscience, natural language processing, and related research fields.

Many biological characteristics arise by interactions between more than one biological organism or unit. Fertilization success in sexually reproducing species represents such an extended phenotype where both mates are required to be fertile for a successful outcome. Consequently, predictive models should account for the joint nature of reproductive performance while offering interpretable estimates for individual mate contributions. Recent advances in genomics and machine learning (ML) provide standardized, high-dimensional genetic information on one hand and computational tools capable of modeling complex biological systems on the other. Here, we construct and evaluate two-tower (TT) machine learning architectures for genomic prediction of binary reproductive outcomes and recovery of sex-specific fertility liabilities. Simulated datasets, generated under a range of genetic architectures, were utilized to compare multilayer perceptron (TT-MLP), convolutional neural network (TT-CNN), and L1-regularized linear (TT-LASSO) two-tower models. Simulation scenarios varied sex-specific heritabilities, genetic correlations, infertility prevalence, mating structure, and sex-specific infertility rates. Models were evaluated with regard to their ability to predict reproductive success at pair level and also recover true underlying genetic values for male and female fertility. Prediction accuracy increased with the underlying heritable component as expected, while sex-specific tower-scores successfully recovered latent fertility liabilities despite models being trained only on observed joint outcomes. TT-LASSO achieved the highest overall classification performance, whereas TT-MLP provided more balanced and consistent recovery of sex-specific genetic values across scenarios. An additional simulation, incorporating genotype-dependent mate compatibility demonstrated advantages of fully-connected neural networks for capturing non-additive interactions. These results indicate that two-tower frameworks provide a powerful approach for modeling reproductive traits, enabling simultaneous prediction of aggregate reproductive outcomes and sex-specific fertility liabilities from genotypic information.

Segregation of speech into syllables is a key step in neural speech processing. It relies on the alignment of neural activity with the rhythmic structure of speech. Two competing hypotheses explain this neural speech tracking, phase-resetting and evoked responses. While phenomenological modelling of these hypotheses has been successful, we still lack understanding of the underlying cortical circuits. To investigate these mechanisms, we evaluate whether a biophysical next-generation neural mass model can reproduce several features of neural speech tracking, using phenomenological models of the competing hypotheses as algorithmic baselines. We investigate the models dynamics with four tests: recreating in-silico an EEG experiment that identified a correlation between tracking strength and phoneme sharpness, computing the Phase Concentration Metric, testing the effect of varying syllabic rates, and evaluating the Inter Event Phase Coherence across phoneme onsets. While all of the models that we study reproduce the sharpness-tuned rhythmic speech tracking, the evoked model requires a pre-processed acoustic edge impulse stimulus. We demonstrate that the neural mass model is performing thresholded phase-resetting triggered by sharp onsets in the continuous speech envelope. This produces cross-frequency nested oscillations that qualitatively match an experimentally-observed dual-peak signature in the Inter Event Phase Coherence. Our results indicate that the biophysical neural mass model provides a mechanistic bridge between generic oscillatory dynamics in cortical populations and the cognitive computations of speech tracking. Indeed, the non-linear dynamics of the neural mass model offer an explanation for how peak-rate event representations in auditory cortex activity arise in response to continuous acoustic input. Significance StatementSyllable segregation is crucial but challenging as natural speech lacks clear boundaries, yet humans perform this computation effortlessly. Speech aligns neural activity to syllabic rhythms, predicting syllable timing, but the underlying cortical mechanisms remain unknown. Relating this macroscopic behaviour to neurobiology is challenging; however, next-generation neural mass models promise to resolve this. We demonstrate that these models reproduce sharpness-tuned tracking and acoustic edge extraction. Dynamical analyses indicate this occurs through thresholded phase-resetting to phoneme onsets, triggering cross-frequency nested oscillations. Our results both advance biophysical understanding of syllable segregation and validate the models capacity for simulating macroscopic neural activity. These models offer a bridge between the neurobiology of the auditory cortex and speech processing dynamics that phenomenological models cannot provide.

Reliable benchmark datasets are critical for advancing EEG-based emotion recognition. The Finer-grained Affective Computing EEG Dataset (FACED) is the largest publicly available EEG emotion dataset (123 subjects, nine emotion categories) and a widely adopted benchmark. We demonstrate that both intra-subject and cross-subject classification on FACED primarily reflects stimulus identity rather than emotion. Using a linear classifier (LinearSVC) and a deep learning model (CLISA), we show that (1) classification performance is comparable for trials where subjects reported feeling versus not feeling the assigned emotion; (2) accuracy drops when stimulus-assigned labels are replaced with individual self-reports; and (3) accuracy increases when reducing to one video per emotion despite discarding two-thirds of the data. These results reflect three design choices in FACED: few stimuli per category, stimulus-assigned labels, and within-video temporal splits for cross-validation. Together, these make the dataset susceptible to temporal autocorrelation and stimulus-identity confounds. To guide future work, we propose five recommendations -- spanning stimulus diversity, temporal independence, and label validation -- for emotion-decoding study designs that mitigate these confounds.

Evaluating the potential applications of a medicine is a fundamental challenge in drug development. There is a lack of standardized, decision-oriented benchmarks that test whether computational models can generalize therapeutic hypotheses across diseases in ways that reflect real-world pharmaceutical investment decision making. To address this gap, we introduce two complementary resources: the Indication Expansion Investment Decision Network (IxIDN) and the Orphanet Rare Disease Ontology Negative-network (ORDON). IxIDN is a clinical-trial-derived positive benchmark constructed by projecting drug-disease associations from pharmaceutical clinical trials into a disease-disease network; each edge connects disease pairs that have entered clinical trials for the same drug, thereby capturing cases when concrete indication-expansion decisions have been made. The current release contains 574 rare diseases and 5,336 edges. In contrast, ORDON serves as a stringent, biology-aware negative benchmark derived from the authoritative Orphanet Rare Disease Ontology. It identifies maximally distant disease pairs according to curated hierarchical structure and genetics-linked inheritance patterns, providing 793 rare diseases and 5,000 edges that represent high-separation negative candidates across therapeutic areas. Together, IxIDN and ORDON enable rigorous cross-evidence generalization from clinical trials to disease ontology, testing for Disease-Disease Association Learning (DDAL), a core task for mechanism-centered drug repurposing and indication expansion. All data are publicly available with detailed metadata, enabling reproducible evaluation of models on transparent, decision-relevant benchmarks.

Artificial Intelligence (AI) frameworks for automating scientific research have shown strong performance on benchmarks, but their capacity to routinely reproduce results from multiple real-life published studies remains largely untested. We evaluated five advanced AI research frameworks (Kosmos, K-Dense, ToolUniverse, BioAgents from bio.xyz, and the AI Scientist-v2 from Sakana AI) on three real-life tasks (including two recently published papers) spanning uncertainty quantification for molecular property predictions, machine learning on Therapeutic Data Commons benchmarks, and agent-based modeling. AI frameworks demonstrated genuine strengths: generating original hypotheses, competently executing routine data acquisition and coding tasks, providing statistical measures of confidence often absent from the original papers, and producing well-formatted final reports. At the same time, our experiments revealed that real-world scientific tasks remain considerably harder than current benchmarks suggest. No AI framework matched the scope or depth of the original studies, results varied across multiple runs of the same framework with the same prompt, and we documented cases of severe hallucinations in final reports, gaps in literature coverage, and overconfident conclusions. Verification of AI outputs required substantial domain expertise. While these three tasks are only partially representative of the broader scientific landscape, they offer a starting point for developing a more rigorous methodology for evaluation of AI performance than what is currently practiced. We conclude that AI frameworks are already valuable for prototyping research directions and stress-testing completed studies, and some of the limitations documented here appear largely tractable through infrastructure improvements and continued development.

Coexpression of transcripts presents the most common means of computational inference of transcription factor regulation, and is often combined with other data types to infer regulatory networks. With the growing popularity of single-cell approaches, there are questions about how best to extract coexpression information from the data. Recently we reported a simulation study that explored the differences among coexpression performed at different levels: across single cells (xCell, per cell type), across subjects from pseudobulked single-cell data (xSubject, per cell type), or across subjects using bulk tissue samples (xBulk). Here we test predictions made by those models using real data. We consider both preservation (consistency of coexpression findings across different levels of analysis of the same data) and replicability across independent studies, as well as biological interpretability. We find that preservation across levels is limited, indicating the choice of analysis level will affect outcomes. We show that xCell coexpression is more replicable across studies compared to xSubject. xBulk coexpression is dominated by patterns driven by variability in cellular composition and fails to capture much coexpression that is reliably detected at finer resolutions. While all modes of analysis exhibit some enrichment for known regulatory relationships, it was highest with the xCell mode. Finally, we present a case study of the effect of analysis modes on a schizophrenia-associated pattern, reinforcing the importance of analytic choices in the interpretation and replicability of coexpression analyses. Together with our modeling study, this work emphasizes the importance of understanding sources of expression covariation as they relate to the goals of the analysis, and recommend single-cell-based data with biological replicates should be the focus of attempts to infer dynamic regulatory interactions that are more likely to be replicable by others.

Cancer-associated signaling pathways often exhibit abnormal activation under simultaneous dysregulation of multiple molecular components. This study presents a probabilistic temporal Dynamic Bayesian Network (DBN)-based framework for analyzing multi-fault behaviour and intervention response in Growth Factor (GF) and Mitogen-Activated Protein Kinase (MAPK) signaling pathways. Unlike deterministic Boolean propagation, the proposed model represents each pathway component through an activation probability and propagates these probabilities over discrete time steps using soft-logic update rules. One-, two-, three-, and four-fault scenarios were systematically evaluated under a common lowest-burden input vector. The resulting output probabilities were summarized using an encoded pathway-burden score, and known-drug combinations were ranked using efficiency scores relative to no-intervention baselines. Pareto analysis was further used to balance intervention efficiency against drug-vector burden, while a custom dual-target search was performed to identify computational intervention hypotheses beyond predefined drug targets. Results showed that encoded burden increased with fault order in both pathways, with MAPK producing a higher baseline burden than GF. Among known-drug vectors, U0126+LY294002+Temsirolimus consistently emerged as the strongest low-burden candidate, achieving efficiency close to the maximum six-drug vector. Custom dual-target analysis identified ERK1/2+RPS6KB1 in GF and Raf+MEK1 in MAPK as high-impact computational target pairs. Runtime benchmarking showed that batched vectorized NumPy execution substantially improved scalability for higher-order fault simulations. Overall, the framework provides an interpretable and scalable approach for probabilistic pathway-level fault analysis and intervention prioritization.

Endometriosis, a common inflammatory gynecological disorder affecting up to 10% of women worldwide, is characterized by the presence of endometrium-like tissue outside the uterus. Current diagnostic methods, such as ultrasound and MRI, effectively detect ovarian and deep endometriosis but fail to detect more common peritoneal type. Diagnosing peritoneal endometriosis currently necessitates invasive laparoscopy and histological confirmation. Despite numerous efforts, no new reliable biomarkers have successfully transitioned into routine clinical use. This study aimed to investigate the use of targeted metabolomics to discover metabolite ratios capable of identifying endometriosis in plasma samples. We analyzed a discovery population of 235 patients and a validation population of 278 patients. All cases and controls in both populations were diagnosed by laparoscopy. Control subjects included individuals presenting with symptoms such as pain, dysmenorrhea, infertility, or other benign conditions, but who had no laparoscopic evidence of endometriosis. Using generalized linear models (GLMs) and machine learning, the study identified specific metabolite ratios as potential biomarkers that can distinguish different types of endometriosis and enable mass spectrometry-based diagnostics for peritoneal endometriosis. The best-validated GLM, derived from the concentration ratios of amino acids, acylcarnitines, sphingomyelins, and phosphatidylcholines, consisted of Thr/SM(OH) C22:2 + PC aa C40:5/SFA_PC + lysoPC a C16:0/SM(OH) C16:1. This model yielded an AUC of 0.82 (95% CI 0.619-0.891, with 76% sensitivity and 81% specificity) for peritoneal endometriosis. This innovative approach offers a robust diagnostic model, addressing an unmet medical need by facilitating earlier detection of peritoneal endometriosis and improving overall clinical management.

Optimizing dose and schedule remains a central challenge in oncology drug development, particularly for immunotherapies where fixed dosing regimens often fail to account for patient specific heterogeneity in tumor-immune dynamics. Here, we present a hybrid quantitative systems pharmacology-reinforcement learning-Monte Carlo Tree Search (QSP-RL-MCTS) framework for personalized immunotherapy dosing that formulates dose selection as a sequential decision-making problem. The approach integrates a mechanistic QSP model of prostate cancer immunotherapy, transcriptomics informed virtual patient populations and data driven AI system comprising reinforcement learning and Monte Carlo tree search. Reinforcement learning is used to learn adaptive generalized dosing policies that optimize treatment outcomes across the population, while Monte Carlo Tree Search provides forward-looking evaluation of RL predicted dosing trajectories to refine patient-specific decisions. On benchmarking against fixed dosing regimens of ipilimumab, the remission rate of the proposed model (95.2%) was comparable to the highest fixed dosing regimen of 10 mg/kg per dose while the median total dose (72 mg/kg) of the proposed model designed regimen was comparable to the lowest fixed dosing regimen of 3 mg/kg per dose. The model is generalizable across different dosing protocols and can be extended to predict optimal dose under different therapeutic scenarios. Analysis of the learned dosing trajectories enables stratification of patients into distinct response groups and identifies drug activity rate as the dominant determinant of long-term treatment outcome. These results demonstrate how mechanistically guided artificial intelligence can transform population-level dose optimization into patient-specific, biologically interpretable treatment strategies for precision immuno-oncology.

Agentic artificial intelligence (AI) systems increasingly claim to automate scientific research, yet independent evaluations report persistent gaps between those claims and demonstrated capability. We tested frontier agentic AI systems on three practical problems: prediction of treatment non-response in immune-mediated inflammatory diseases, optical chemical structure recognition for literature mining, and prediction of drug-design-related properties from small datasets. Each problem was first assigned to autonomous frameworks and then reattempted as human-led, AI-assisted work. Autonomous runs failed in most cases, while human-led work produced reusable resources and modest but defensible performance, including new evidence for possible mechanisms of treatment resistance and a more practical benchmark for mining chemical structures from scientific papers. Property prediction was the single task on which one autonomous AI framework matched the human expert. We conclude that current frameworks can carry out engineering and analysis once a human expert leads the project, but cannot yet engineer a novel solution without oversight. The use of AI on real-life scientific problems remains an art rather than a routine technology.

Background/ObjectivesMetabolism-level interpretation of metabolomics datasets requires aggregation analyses across metabolites. One highlyused aggregation analysis is pathway enrichment analysis (PEA), which involves detecting pathways enriched with metabolites that are differential between experimental groups. Annotating metabolites with pathway associations is a prerequisite for PEA. While several knowledgebases define pathways and include metabolite-pathway annotations, these definitions are often partially or even grossly incomplete due to limitations in current metabolic knowledge and its curation, which greatly limits the effectiveness of PEA. MethodsIn this work, we used a novel multitask classification, graph convolutional-like neural network to generate high-quality metabolite-pathway annotations for pathways defined across KEGG, MetaCyc, and Reactome. We then included these predicted metabolite-pathway annotations when performing PEA on 990 datasets deposited in Metabolomics Workbench. ResultsWe demonstrate an 8-fold increase in the median number of enriched pathways detected across these datasets compared to using only knowledgebase-derived annotations. ConclusionsThe significant increase in enriched pathways substantially improves the biological and biomedical interpretability of metabolomics datasets.

Predicting gene essentiality across cellular contexts is a central challenge in computational biology, with implications for identifying cancer vulnerabilities. Graph neural networks (GNNs) integrate molecular interaction networks with gene-level features, but it remains unclear whether their performance gains arise from biologically meaningful connectivity or generic graph structure. Here, we systematically evaluate the role of network information in gene essentiality prediction using 2,741 genes across three tissues. We compare GNNs to feature-only baselines, including multilayer perceptron (MLP) and random forest (RF) methods, under a strict gene-level 5-fold cross-validation scheme to prevent information leakage. To isolate the role of network information, we assess models on the STRING protein-protein interaction network, a degree-preserving shuffled network, and a fully randomized network, with and without network-derived features. GNNs outperform feature-only models, reducing mean squared error and improving Matthews correlation coefficient across all tissues. However, these gains depend critically on biologically structured connectivity: performance degrades substantially under randomized topology and is not preserved by degree-constrained rewiring. Network features are largely redundant when using biologically meaningful graphs, as their information is recovered through message passing, but become important when topology is uninformative. Per-gene analyses reveal uniformly low correlations across models, highlighting intrinsic limits imposed by data variability. Graph Transformer models incorporating global attention do not outperform standard GNNs, indicating that predictive signals are predominantly local. Together, these results show that predictive gains arise from biologically structured connectivity rather than generic graph topology.

Mathematical models of enzyme cycles form the basis of quantifying key features of metabolism and membrane transport. These models are often integrated into more comprehensive models such as whole-cell models to understand emergent behaviours between interacting components. However, it is currently computationally infeasible to simulate the full dynamical behaviour of every enzyme at a network scale. Model reduction is frequently used to improve computational efficiency, but in general, these approaches do not preserve physical and thermodynamic consistency. Here, we outline a general method for simplifying enzyme kinetics models while retaining mass, charge and energy balance. We base our approach on the bond graph, which is a general methodology for modelling biological systems from fundamental physical laws. This approach ensures that key physical constraints are enforced in every model, regardless of their complexity. Our thermodynamic model reduction framework is readily extended to electrogenic transporters through the coupling of chemical and electrical processes. Through the application of our approach to both hypothetical enzyme cycles and real data from the Na+/K+ ATPase, we show that it can rapidly screen for plausible network structures in circumstances where enzyme catalytic mechanisms may not be fully characterised, facilitating biological discovery and drug development.

Biological databases store curated knowledge that researchers traditionally access through web interfaces or APIs. To move beyond casual browsing requires domain-specific knowledge and expertise to frame the queries necessary to explore this data. This generates a barrier for new users in scientific fields undergoing paradigm shifts. Exposing these databases to large language models (LLMs) via the Model Context Protocol (MCP) enables natural-language access, a potential accessibility solution. We implement this for Virtual Fly Brain (VFB), an expert-curated and ontology-backed knowledgebase of Drosophila neuroscience, providing the precision needed to make recently-integrated connectomes accessible. Benchmarked on 30 neuroscience tasks against a bare LLM and a web-search-assisted LLM, the VFB-MCP-equipped LLM produces precise, verifiable and appropriately quantified answers on 25/30 tasks vs 14/30 for web and 2/30 for bare (Wilcoxon p<0.01, Holm-corrected, all pairwise comparisons). The MCP advantage is largest for tasks where data quantification is required (89% vs 11% web). This work establishes MCP over ontology-backed knowledge graphs as an effective method to improve LLM response quality for neuroscience and connectomics data.

Metabolites from Plantago lanceolata (plantain) biomass have been linked with biological nitrification inhibition (BNI) in soil. After grazing, leaf metabolite chemistry is altered via digestion, and a suite of secondary metabolites are then delivered onto soil via dung and urine. The purpose of this study was to establish if urine from sheep grazed on plantain had BNI activity when added to pasture soil, and to identify the metabolite profile(s) that most likely contribute to the BNI effects observed. Groups of sheep (n=5) were grazed on one of nine different plantain cultivars in autumn and spring with analysis of leaf material, urine, soil incubation and BNI bioassay data used to identify potential metabolite candidates implicated with BNI. The urinary nitrogen and metabolite composition of sheep fed plantain varied significantly between cultivars and season. After 28 days of incubation, all soil microcosms treated with plantain-derived urine had up to 35% less nitrate than comparative ryegrass urine controls in both seasons, except one in autumn. The key phytochemistry associated with lower soil nitrate concentrations was phenylethanoid and iridoid glycosides resulting in a higher output of glucuronidated, methylated and sulfated secondary metabolites in the urine. Among 19 secondary metabolites identified in the urine, hydroxytyrosol-related metabolites as well as catechol glucuronide, 2-methoxyphenyl sulfate and guaiacol-{beta}-D-glucuronide appear to be the most likely target compounds with respect to the BNI effects observed. Variation in metabolites from different plantain cultivars affected the ratio of metabolite derivatives in urine, which ultimately affected soil nitrification rates. Cultivar phytochemistry is therefore an important consideration with respect to BNI under urine patches. HighlightsO_LISheep grazing different plantain cultivars had different urine compositions C_LIO_LIUrines elicited biological nitrification inhibition (BNI) in soil and in vitro C_LIO_LIDifferent BNI response was related to differential expression of urine metabolites C_LIO_LIKey urine metabolites associated with BNI are derived from glycosidic compounds C_LI